##########################################################################################

library(dplyr)
library(data.table)
library(optparse)
library(reshape)

##########################################################################################

option_list <- list(
    make_option(c("--sig_file"), type = "character") ,
    make_option(c("--reprot_file"), type = "character") ,
    make_option(c("--mutRate_file"), type = "character") ,
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    info_file <- "~/20230331_colrNew/public_data/liver/config/tumor_normal.class.list"
    maf_path <- "~/20230331_colrNew/public_data/liver/results/maf"
    out_path <- "~/20230331_colrNew/public_data/liver/results/mutation_QC"

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sig_file <- opt$sig_file
reprot_file <- opt$reprot_file
mutRate_file <- opt$mutRate_file
out_path <- opt$out_path

dir.create(out_path , recursive = T)
dir.create( paste0( out_path , "/Input_all_VAF005_read4"  ))

###########################################################################################

info <- fread( info_file ,data.table=F)

###########################################################################################
## ------------------------------------------------ ##
## Quality control on mutations prior to Pyclone
## ------------------------------------------------ ##
# VAF≥0.05 and variant read count>3
# VAF≥0.01 and variant read count>3

for(normal in unique(info$Normal)){

  sinf <- subset(info, Normal == normal)
  sinf$Class_a <- sinf$Class_sub

  # read in maf files
  maf_file <- paste(maf_path , "/" ,sinf$Tumor[1],'_',sinf$Normal[1],'_GGA_Filter_funcotated.maf',sep='')
  gga <- read.csv( maf_file, comment.char = "#", sep='\t')
  
  gga$mutID <- paste(gga$Hugo_Symbol, gga$Variant_Classification, gga$Chromosome, gga$Start_Position, gga$Reference_Allele, gga$Tumor_Seq_Allele2, sep=':')

  print( table(duplicated(gga$mutID)) )
  
  # calculate Depth and VAF
  dat <- gga[,c('mutID','t_alt_count','t_ref_count')]
  dat$Depth <- dat$t_alt_count + dat$t_ref_count
  dat$VAF <- dat$t_alt_count / dat$Depth
  
  dat$qc1 <- ifelse( dat$VAF >= 0.01 & dat$t_alt_count > 3, 1, 0)
  dat$qc2 <- ifelse( dat$VAF >= 0.05 & dat$t_alt_count > 3, 1, 0)
  
  colnames(dat)[2:ncol(dat)] <- paste(colnames(dat)[2:ncol(dat)],sinf$Class_a[1],sep='.')

  data <- dat
  rm(gga, dat)  

  if(nrow(sinf) >= 2){ 
    for(i in 2:nrow(sinf)){

      # read in maf files
      maf_file <- paste(maf_path , "/" ,sinf$Tumor[1],'_',sinf$Normal[1],'_GGA_Filter_funcotated.maf',sep='')
      gga <- read.csv(maf_file, comment.char = "#", sep='\t')
    
      gga$mutID <- paste(gga$Hugo_Symbol, gga$Variant_Classification, gga$Chromosome, gga$Start_Position, gga$Reference_Allele, gga$Tumor_Seq_Allele2, sep=':')

      print( table(duplicated(gga$mutID)) )
    
      # calculate Depth and VAF
      dat <- gga[,c('mutID','t_alt_count','t_ref_count')]
      dat$Depth <- dat$t_alt_count + dat$t_ref_count
      dat$VAF <- dat$t_alt_count / dat$Depth
    
      dat$qc1 <- ifelse( dat$VAF >= 0.01 & dat$t_alt_count > 3, 1, 0)
      dat$qc2 <- ifelse( dat$VAF >= 0.05 & dat$t_alt_count > 3, 1, 0)
    
      colnames(dat)[2:ncol(dat)] <- paste(colnames(dat)[2:ncol(dat)],sinf$Class_a[i],sep='.')

      data <- merge(data, dat, by='mutID', all=T)
      rm(gga,  dat, i)  
    }
  }

  out_file <- paste( out_path , "/" , normal,'_prePyclone_QC_note.maf',sep='')
  write.table(data, file = out_file, row.names=F, quote=F, sep='\t')

  rm(data)
  rm(normal)

}


###########################################################################################
# 质控标准v2 ----------------------------------------------#
# (1) 若sSNV不满足 VAF≥0.01 and variant read count>3，则将其 t_alt_count 转为0
# (2) 若无样本满足 VAF≥0.05 and variant read count>3，则去除该突变
for(normal in unique(info$Normal)){

  pri <- read.delim( paste( out_path , "/" , normal,'_prePyclone_QC_note.maf',sep='') )
  
  pri$qc1_sum <- apply(pri[, colnames(pri)[grep('qc1.',colnames(pri))] ,drop=F], 1, function(x){return(sum(x,na.rm=T))})
  table(pri$qc1_sum)
  
  pri$qc2_sum <- apply(pri[, colnames(pri)[grep('qc2.',colnames(pri))] ,drop=F], 1, function(x){return(sum(x,na.rm=T))})
  table(pri$qc2_sum)
  dat <- subset(pri, qc2_sum > 0)
  rm(pri)
  
  # 若sSNV不满足 VAF≥0.01 and variant read count>3，则将其 t_alt_count 转为0
  sam <- info[info$Normal==normal,'Class_a']
  sam <- gsub('-','_',sam)
  for(j in sam){
    altname <- paste('t_alt_count',j,sep='.')
	qcname <- paste('qc1',j,sep='.')
	dat[,paste(altname,'v1',sep='.')] <- ifelse(dat[,qcname] == 0, 0, dat[,altname])
	rm(altname, qcname, j)
  }
  rm(sam)
  
  # 注释SNV信息，准备Pyclone的输入文件
  case_info <- info[info$Normal==normal,]
  for(i in 1:nrow(case_info)){
  	chat <- read.delim(paste('/public/user/wyz2014/colr_mutiple/titan/chat/',case_info$Tumor[i],'_',case_info$Normal[i],'_GGA_CHAT.txt',sep=''))
  	
  	chat$mutID <- paste(chat$Hugo_Symbol, chat$Variant_Classification, chat$Chr, chat$Start_Position, chat$REF, chat$ALT, sep=':')
  	
  	samaf <- dat[,c('mutID', colnames(dat)[grep( gsub('-','_',case_info$Class_a[i]),colnames(dat))] )]
  	colnames(samaf) <- gsub( paste('.',gsub('-','_',case_info$Class_a[i]),sep=''),'',colnames(samaf))
  	
  	mafcn <- merge(samaf,chat[,c('mutID', 'minor_cn', 'total_cn')],by='mutID',all.x=T)
  	mafcn <- subset(mafcn, !is.na(VAF))
  	mafcn$mutation_id <- mafcn$mutID
  	mafcn$ref_counts <- mafcn$t_ref_count
  	mafcn$var_counts <- mafcn$t_alt_count.v1	
  	mafcn$normal_cn <- 2
  	mafcn$major_cn <- mafcn$total_cn - mafcn$minor_cn
  	
  	infile <- mafcn[,c('mutation_id','ref_counts','var_counts','normal_cn','minor_cn','major_cn')]
  	infile$minor_cn <- ifelse(is.na(infile$minor_cn),1,infile$minor_cn)
  	infile$major_cn <- ifelse(is.na(infile$major_cn),1,infile$major_cn)
  	infile <- subset(infile, major_cn>0)
  	
  	write.table(infile,file=paste('/public/user/wyz2014/CRC_metastasis/5_clone_evolution/Input_all_VAF005_read4/',normal,'/',case_info$Class_a[i],'.tsv',sep=''),row.names=F,quote=F,sep='\t') # 请核查路径
  	
  	rm(chat, samaf, mafcn,infile,i)
  }
  rm(case_info, normal,dat)
}


